G x E expression

The question is how much G x E we see in the RNAseq data.

To get started on this I ran the file “ProcessCount.R”. For this the data was TMM normalized, voom transformed, and then a gtXenvironment model was fit in limma.

library(limma)
library(tidyverse)

── Attaching packages ─────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 2.2.1     ✔ purrr   0.2.4
✔ tibble  1.3.4     ✔ dplyr   0.7.4
✔ tidyr   0.7.2     ✔ stringr 1.2.0
✔ readr   1.1.1     ✔ forcats 0.2.0
── Conflicts ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(stringr)
library(modelr)
library(broom)

Attaching package: ‘broom’

The following object is masked from ‘package:modelr’:

    bootstrap

library(multidplyr)

Get the limma objects

load("voom.int.Rdata")

How many genes show GxE?

generate vector of interaction coefficients and then topTable

colnames(fit.int$coefficients)
int.coef.n <- colnames(fit.int$coefficients) %>% str_which(":") 
topTable(fit.int,coef = int.coef.n)

Not working

let’s just do an ANOVA for each gene

counts <- as.tibble(counts.voom.gr$E) %>% 
  mutate(geneIndex=1:n())
counts[1:10,1:10]

counts.long <- counts %>% gather(key="name",value="expression", -geneIndex) %>%
  mutate(gt=factor(str_extract(name,"RIL_[0-9]*")), 
         env=factor(str_extract(name,"CR|UN")), 
         rep = factor(str_extract(name,"Rep[0-9]"))) %>%
  select(-name) %>%
  group_by(geneIndex) %>% nest()
counts.long

counts.long$data[[1]]

Now fit an anova and do a VCA for each gene in parallel

fitAnova <- function(df) aov(expression ~ gt*env,data=df)
getPvals <- function(fit) {
  pvals <- summary(fit)[[1]]$`Pr(>F)`
  names(pvals) <- dimnames(summary(fit)[[1]])[[1]] %>% str_trim()
  data.frame(t(pvals))
}
fitVCA <- function(df) {
  tmpvca <- anovaVCA(expression ~ gt*env,Data=as.data.frame(df))
  tmpvca$aov.tab[,"%Total"]
  }
  
cl <- create_cluster(cores=4)

Initialising 4 core cluster.

counts.long.test <- partition(counts.long[1:100,],geneIndex,cluster = cl)

group_indices_.grouped_df ignores extra arguments

cluster_library(cl,c("purrr","broom","VCA","stringr"))
cluster_copy(cl,fitAnova)
cluster_copy(cl,getPvals)
cluster_copy(cl,fitVCA)
system.time(
  counts.long.test <- counts.long.test %>% mutate(aov = map(data,fitAnova), 
                                                  glance = map(aov,glance),
                                                  pvals = map(aov,getPvals),
                                                  varcomp = map(data,fitVCA))
)

   user  system elapsed 
  0.039   0.041  17.603